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ABSTRACT 

Phase  transitions  occur  in  a  variety  of  physical  phenomena  ranging  from  condensation  of 
gases  to  solidification  of  liquids.  It  is  of  paramount  importance  to  provide  an  adequate 
modelling  of  these  phenomena  in  a  thermodynamically  consistent  way.  This  paper  ad¬ 
dresses  an  introduction  into  the  mathematical  modelling  of  phase  transitions  in  fluids  from 
the  perspective  of  consistently  employing  a  modified  Legendre  transform  of  entropy  con¬ 
sidered  as  a  given  function  of  internal  energy  and  volume.  An  explicit  conservative  scheme 
for  incompressible  phase  transition  and  an  implicit  non-conservative  scheme  based  on  the 
fictitious-capacity  and  lubrication- theory  approximations  are  implemented  in  MAT  LAB©. 
Illustrative  numerical  simulations  are  conducted,  and  some  results  are  verified  against,  an 
exact  benchmark  solution. 
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Executive  Summary 

Phase  transitions  occur  in  a  variety  of  physical  phenomena  ranging  from  condensation  of 
gases  to  solidification  of  liquids  at  certain  conditions.  It  is  of  paramount  importance  to 
provide  an  adequate  modelling  of  these  phenomena  in  a  thermodynamically  consistent  way. 
This  paper  addresses  an  introduction  into  the  mathematical  modelling  of  phase  transitions 
in  fluids  from  the  perspective  of  consistently  employing  a  modified  Legendre  transform  of 
entropy  considered  as  a  given  function  of  internal  energy  and  volume  (primary  variables) . 
The  Legendre  transform  depends  naturally  on  the  absolute  temperature  and  pressure 
(potential  variables),  and  coincides  with  the  Gibbs  free  energy.  When  entropy  ceases  to  be 
a  strictly  concave  function  of  the  primary  variables,  the  derivatives  of  the  Gibbs  free  energy 
become  discontinuous  with  respect  of  the  potential  variables  at  some  phase-transition 
line.  These  derivatives  define  the  primary  variables  at  equilibrium,  and  their  jumps  are 
the  latent  energy  and  volume  of  the  phase  transition.  The  second-order  derivatives  of 
the  Gibbs  free  energy,  defining  heat  capacity,  compressibility  and  expansivity,  exhibit 
delta- function  singularities  at  the  phase-transition  line.  The  main  feature  of  a  Legendre 
transform  is  concavity  (or  convexity)  from  which  a  bunch  of  important  consequences  are 
readily  obtained. 

An  explicit  conservative  scheme  for  incompressible  phase  transition  and  an  implicit  non¬ 
conservative  scheme  based  on  the  fictitious-capacity  and  lubrication-theory  approxima¬ 
tions  are  implemented  in  MATLAB®.  Illustrative  numerical  simulations  are  conducted, 
and  some  results  are  verified  against  an  exact  benchmark  solution. 
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1  Introduction 

Phase  transitions  occur  in  a  variety  of  physical  phenomena  ranging  from  the  condensation 
of  gases  to  solidification  of  liquids  at  certain  conditions.  It  is  of  paramount  importance  to 
provide  an  adequate  mathematical  modelling  of  these  phenomena  in  a  thermodynamically 
consistent  way.  The  main  objective  of  this  paper  is  to  present  an  introduction  into  the 
thermodynamics  of  phase  transitions  using  the  useful  mathematical  apparatus  of  Legendre 
transform  [Zia,  Redish  Sz  McKay  2009],  formulate  mathematical  models  starting  from  the 
most  general  one,  and  describe  numerical  schemes.  Illustrative  numerical  simulations  are 
conducted,  and  some  schemes  are  verified  against  an  exact  solution. 

The  Legendre  transform  is  commonly  employed  in  Classical  Mechanics  [Arnold  1989]  and 
Thermodynamics  [Callen  1985],  however  its  power  is  not  always  utilized  in  full  as  the 
traditional  approach  of  teaching  Physics  is  historically  dominated.  We  will  show  here  how 
many  important  consequences  can  be  derived  economically  in  a  straightforward  manner 
from  the  proper  definition  of  Gibbs  free  energy  employing  the  Legendre  transform  whose 
use  is  inextricably  linked  with  the  Maximum  Entropy  Principle. 

For  the  sake  of  simplicity,  we  confine  ourselves  to  phase  transitions  in  fluids  composed  of 
a  single  species,  which  comprise  liquids  and  their  vapours  (gases),  and  can  be  modelled 
by  a  two-parametric  thermodynamic  system.  The  latter  means  that  only  two  independent 
variables  define  the  thermodynamic  state  of  the  system,  such  as  energy  and  volume,  the 
primary  variables,  or  temperature  and  pressure,  the  potential  variables.  An  extension  of 
this  approach  to  fluids  composed  of  several  species  is  straightforward. 

Though  both  types  of  variables  are  generally  equivalent  for  the  description  of  a  physical 
phenomenon,  the  transformation  from  the  primary  to  potential  variables  is  not  always 
invertible.  This  particularly  happens  when  phase  transitions  take  place.  Indeed,  a  small 
variation  in  temperature  and  pressure  across  the  phase-transition  line  in  the  phase  diagram 
results  in  large  changes  in  energy  and  volume.  These  changes  are,  respectively,  the  latent 
energy  and  latent  volume  of  the  thermodynamic  system  undergoing  a  phase  transition. 
For  an  incompressible  phase  transition,  the  latent  volume  vanishes  and  the  latent  energy 
becomes  equal  to  the  latent  heat. 

The  typical  phase  diagram  in  the  (T,  p)-plane,  where  T  denotes  the  absolute  temperature 
and  p  pressure,  is  shown  in  Figure  1.  Two  special  states,  called  respectively  the  critical 
point  and  triple  point ,  are  depicted.  The  critical  point  is  the  point  beyond  which  there  is 
no  distinction  between  liquid  and  its  vapour.  In  other  words,  in  the  supercritical  region 
there  is  a  single  fluid  phase  with  continuous  dependence  of  energy  and  volume  on  pressure 
and  temperature.  The  triple  point  is  the  unique  pressure  and  temperature  for  which  three 
phases  of  a  substance  coexist.  The  sublimation,  fusion  and  vaporisation  lines  branching 
from  the  triple  point  are  shown.  Note  that  additional  phase  transitions  can  take  place  in 
the  solid  state,  which  cannot  be  described  by  a  two-parametric  thermodynamic  system. 

In  many  practical  applications  the  solid-liquid  phase  transitions  can  be  assumed  incom¬ 
pressible,  and  therefore  the  corresponding  thermodynamic  system  becomes  one-parametric. 
Indeed,  the  densities  of  the  two  phases  do  not  change  dramatically,  and  this  is  reflected 
by  the  fact  that  the  fusion  line  is  nearly  vertical  for  most  materials.  This  simplifying 
assumption  allows  one  to  model  incompressible  phase  transitions  by  the  solution  of  the 
Stefan  problem  [Meirmanov  1992]  which  is  decoupled  from  the  fluid  motion.  The  latter 
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Figure  1:  Schematic  of  a  phase  diagram 


means  that  the  velocity  field  can  be  set  zero,  and  only  the  balance  of  energy  constitutes 
the  Stefan  model.  In  contrast  to  this  the  liquid-gas  phase  transitions  occurring  at  the 
vapourisation  line  are  more  complex  in  nature  as  the  corresponding  thermodynamic  sys¬ 
tem  is  genuinely  two-parametric,  and  there  is  no  way  to  decouple  fluid  motion  [Kutateladze 
&  Borishanskii  1966,  Kandlikar  1999].  The  model  becomes  more  complicated  if  addition¬ 
ally  the  effect  of  surface  tension  is  taken  into  account  [Antanovskii  1996]. 


2  Elements  of  thermodynamics 

For  ready  reference,  we  recall  the  basic  notions  and  principles  of  thermodynamics  [Fermi 
1937,  Reif  1965,  Kittel  1969].  The  main  attention  is  paid  to  the  rigorous  definition  of  all 
the  necessary  variables  starting  with  entropy  dependence  on  macroscopic  primary  variables 
and  employing  a  Legendre  transform.  Then  we  discuss  the  thermodynamics  of  phase 
transitions  and  provide  an  illuminating  example  using  the  van  der  Waals  equation  of 
state. 

2.1  General  principles 

Fluids  composed  of  a  single  species,  such  as  pure  water  in  either  liquid  or  vapour  phase, 
can  be  described  by  a  two-parametric  thermodynamic  system.  This  means  that  two  inde¬ 
pendent  variables,  the  internal  energy  E  and  volume  V .  define  the  state  of  a  substance  of 
some  fixed  mass  M  or,  equivalently,  of  some  fixed  number  of  particles  N.  These  parame¬ 
ters  are  reliably  controlled  in  an  experiment,  and  a  system  with  given  E  and  V  is  called 
an  isolated  thermodynamic  system.  In  other  words  the  isolated  thermodynamic  system  is 
not  allowed  to  exchange  energy  with  the  environment  or  do  any  work. 

The  primary  variables  ( E ,  V)  define  the  macrostate  of  the  thermodynamic  system.  How¬ 
ever,  these  variables  do  not  determine  the  system’s  microstates.  It  was  the  greatest  idea 
by  Rudolf  Clausius,  one  of  the  founders  of  Thermodynamics ,  to  introduce  entropy  S  and 
make  it  responsible  for  determining  the  most  probable  microstates  compatible  with  the 
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given  macrostate  through  the  celebrated  Maximum  Entropy  Principle.  A  particular  ma¬ 
terial  has  its  unique  dependence  of  entropy  on  the  primary  variables.  In  other  words, 
to  the  first  approximation,  entropy  is  assumed  to  be  the  only  additional  variable  that 
defines  the  microscopic  structure  of  a  physical  material.  In  particular,  the  dependence 
S  =  S(E,V )  is  capable  of  describing  all  the  observable  physical  phenomena  in  fluids  in¬ 
cluding  phase  transitions.  In  this  sense  Thermodynamics  is  a  very  economic  theory  with 
excellent  applications  to  the  real  world. 

In  this  connection  it  is  worthwhile  noting  that  a  confusion  commonly  arises  from  the  naive 
interpretation  of  the  Maximum  Entropy  Principle  by  stating  that  the  entropy  of  an  isolated 
thermodynamic  system  tends  to  a  maximum.  This  statement  obviously  contradicts  the 
dependence  S  =  S(E,V )  which  tells  us  that  the  entropy  does  not  have  a  freedom  to 
achieve  either  a  maximum  or  a  minimum.  Indeed,  it  is  completely  determined  by  the 
macrostate  through  the  given  dependence  S  =  S(E,V).  This  motivates  the  inclusion  of 
the  rigorous  formulation  of  the  Maximum  Entropy  Principle  below  to  avoid  confusion. 

The  modern  definition  of  entropy  due  to  Ludwig  Boltzmann  is  given  by  the  following 
expression 

S  =  ks  log  W  (1) 

where  ks  is  Boltzmann’s  constant,  W  is  the  number  of  microstates  consistent  with  the 
given  macrostate  (also  called  the  degeneracy  of  the  energy  level),  and  the  logarithm  is 
natural.  This  expression  is  very  useful  in  Statistical  Mechanics ,  the  theory  aimed  at 
deriving  the  macroscopic  properties  of  a  material  from  a  knowledge  of  molecule  interaction, 
that  boils  down  to  calculating  the  degeneracy  W  of  quantum  states  as  a  function  of  E 
and  V .  The  greater  the  degeneracy  W  the  less  information  about  the  system’s  micro¬ 
state  at  particular  instant  is  available,  so  entropy  measures  the  missing  information  about 
the  system.  It  is  postulated  in  Thermodynamics  and  Statistical  Mechanics  that  all  the 
microstates  are  equally  probable. 

The  reason  for  using  the  logarithm  in  the  definition  of  entropy  will  be  clear  from  the 
following  thought  experiment  [Kittel  1969].  Let  us  consider  two  isolated  thermodynamic 
systems  T>a  and  Eb  depicted  in  Figure  2.  The  corresponding  numbers  of  microstates  are 
well  defined  in  terms  of  energies  and  volumes.  Let  T>c  be  the  logical  union  of  the  systems  Va 
and  Eb ,  which  is  also  isolated.  In  this  case  Ec  =  Ea  +  E b,  Vc  =  Va  +  Vj,,  Wc  =  Wa  I'D, 
and  therefore  Sc  =  Sa  +  <S&.  In  other  words,  entropy  behaves  like  an  additive  (extensive) 
variable,  such  as  energy  or  volume. 

Let  us  now  allow  the  two  thermodynamic  systems  to  exchange  energy  and  momentum 
(work)  but  not  mass.  This  means  that  we  are  not  maintaining  the  energy  and  volume 
of  each  open  system,  but  the  total  energy  and  volume  of  the  combined  isolated  system 
Ec  =  Ea  U  Ei).  The  total  energy  and  volume  are  equal  to  the  sums  of  the  initial  energies 
and  volumes  of  the  subsystems  Ea  and  Eb  before  contact  provided  that  the  interfacial 
energy  is  negligibly  small.1  In  other  words, 

Ec  =  E°a  +  E°b,  Vc  =  Va°  +  Vb° ,  (2) 

where  (E°,  V°)  and  ( Eb ,  Vb  )  are  the  initial  energies  and  volumes  of  the  subsystems  Ea 
and  T>b,  respectively. 

lrThis  is  not  the  case  when  the  effect  of  capillarity  is  taken  into  account. 
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Figure  2:  Two  thermodynamic  systems 


The  number  of  microstates  of  the  combined  system  is  equal  to  the  following  sum  of  the 
products 

Wc  ( Ec ,  Vc)  =  Y,  Wa  (Fa,  Va)  Wb  (Eb,  Vb)  (3) 

where  the  summation  is  taken  over  all  the  possible  values  of  energy  and  volume  satisfying 
the  constraints 

Ea  +  Eb  =  Ec  (  =  E°a  +  E°b)  ,  Va  +  Vb  =  Vc  (  =  Va°  +  Vb°)  .  (4) 

Obviously,  the  product  Wa  (E°,V°)  Wb(Eb,Vb)  is  present  in  the  sum  (3)  because  the 
above  constraints  admit  it,  and  therefore  the  total  entropy  satisfies  the  strict  inequality 

Sc  (Ec,  Vc)  >  Sa  (. E°a,  V°)  +  Sb  (E°h,  Vb°)  .  (5) 

This  fact  is  not  surprising  since  two  independent  sets  of  constraints  for  the  thermodynamic 
systems  are  replaced  with  the  single  set  of  constraints  (4)  for  the  combined  system,  and 
therefore  the  number  of  microstates  of  the  combined  system  increases  due  to  extra  freedom. 

The  Maximum  Entropy  Principle  states  that  a  thermodynamic  equilibrium  is  attained  at 
those  pairs  (Ea,  Va)  and  (Eb,  Vb)  which  provide  the  maximum  value  of  some  product  in  (3). 
In  other  words,  we  arrive  at  the  maximisation  problem 

Sa(Ea,Va)  +  Sb(Eb,Vh)->  max  (6) 

under  the  constraints  (4),  which  immediately  gives  the  necessary  conditions  for  thermo¬ 
dynamic  equilibrium2 

dSa(Ea,Va)  =  dSb(Eb,yb)  1  dSg  (Eg,  Va)  =  dSb(Eb,yh )  p 

dE  dE  ~  T  ’  dV  dV  ~  T '  U 

Here  the  definitions  of  the  absolute  temperature  T  and  pressure  p  are  introduced  through 
the  derivatives  of  entropy,  which  are  commonly  called  potential  variables  conjugate  to 
energy  E  and  volume  V.  The  potential  variables  (T,  p)  are  equal  in  each  thermodynamic 
system  in  contact  thus  providing  the  condition  for  thermodynamic  equilibrium. 

2Though  W  is  integer,  we  assume  that  the  resulting  function  S  =  S(E ,  V )  can  be  differentiated. 
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For  given  functions  Sa(E,  V )  and  Sb(E,  V ),  four  unknowns  (Ea,  Va,  Eb,  14)  can  be  in  prin¬ 
ciple  determined  from  the  four  equations  (4)  and  (7).  Alternatively,  each  pair  (Ea,  Va )  and 
(Eb,  Vb)  can  be  determined  in  terms  of  T  and  p  if  these  are  given.  Such  a  situation  arises 
when  one  thermodynamic  system  (T>a)  is  significantly  smaller  than  the  other  one  (74). 
The  larger  system  plays  the  role  of  a  reservoir  of  heat  and  work.  Assuming  Ea  <C  Eb  and 
Va  <  Vb,  the  problem  (6),  equivalently  written  in  the  form 

5a  (Ea,  Va)  +  Sb  (Ec  -  Ea,  Uc  -  Va)  ->•  max ,  (8) 

reduces  to  the  approximate  problem 

Sa  (Ea,  Va)  -  Ea  +TP  Va  ->•  max  .  (9) 

Here  in  the  Taylor  series  expansion  of  the  reservoir  entropy  Sb  (Ec  —  Ea,  Vc  —  Va)  we  retain 
only  the  linear  terms  and  drop  the  constant  term  Sb  (Ec,  Vc)  which  obviously  does  not  affect 
the  solution.  In  this  case  the  absolute  temperature  T  and  pressure  p  are  the  properties 
of  the  reservoir,  which  are  assumed  constant  regardless  of  small  exchange  of  energy  and 
volume  with  the  attached  thermodynamic  system.  Of  course,  the  necessary  conditions  for 
the  solution  of  the  problem  (9)  are  as  follows 

dSa  (Ea,  Va)  _  1  8Sa  (Ea,  Va)  _P 

dE  T  ’  dV  T  ’  1  ] 

which  are  in  agreement  with  (7). 

It  is  worthwhile  emphasizing  that,  for  a  thermodynamic  system  composed  of  a  huge  number 
of  particles,  the  following  equality  holds  with  a  very  high  accuracy 


Sc  (Ec,  Vc)  =  Sa  (Ea,  Va)  +  Sb  (Eb,  Vb)  . 


(11) 


The  approximate  equality  (11)  takes  place  for  thermodynamic  subsystems  in  contact  only 
after  the  equilibrium  is  achieved.  This  is  in  contrast  with  the  strict  inequality  (5). 

This  purely  statistical  observation  allows  one  to  assume  that  entropy  is  an  extensive  ther¬ 
modynamic  variable.  From  this  point  it  makes  sense  to  talk  about  the  specific  energy  e, 
volume  v  and  entropy  s  defined  by  the  formulae 

E  V  S 

e  =  —  ,  v  =  —  ,  s  =  —  . 

M  M  M 


(Alternatively,  we  can  introduce  molar  properties 


e'  = 


NAE 


v  = 


NAV 


N  N 

where  Na  is  Avogadro’s  number.) 


s'  = 


NAS 

N 


(13) 


The  equation  of  state  is  now  reformulated  in  terms  of  the  specific  thermodynamic  prop¬ 
erties:  s  =  s(e,v).  The  subscripts  identifying  thermodynamic  subsystems  are  dropped 
from  this  point.  The  definitions  of  temperature  and  pressure  do  not  change  as  they  are 
intensive  variables 


ds  (e,  v) 
de 


1 

T 


ds  (e,  v) 
dv 


P_ 

T 


(14) 
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The  second  law  of  Thermodynamics  immediately  follows  from  the  above  definitions 

Tds  =  de  +  pdv.  (15) 

This  invariant  form  of  the  second  law  does  not  imply  the  original  set  of  independent 
variables  any  more. 

It  is  worthwhile  stressing  that  the  absolute  temperature  T  is  always  positive.  As  demon¬ 
strated  in  Statistical  Mechanics ,  the  latter  is  due  to  the  fact  that  the  spectrum  of  energy 
of  fluids  does  not  have  an  upper  bound.  In  particular,  the  maximum  entropy  principle  (9) 
is  equivalent  to  the  following 

e  +  pv  —  T  s(e,  v)  — >  min  (16) 

where  the  (specific)  primary  variables  e  and  v  are  free  but  the  potential  variables  T  and  p 
are  fixed.  Of  course,  at  the  point  of  minimum  we  get  (14). 

In  many  textbooks  on  Thermodynamics  the  (specific)  Gibbs  free  energy  is  defined  by  the 
expression 


g  =  e  +  pv  —  T  s  (17) 

which  is  poorly  motivated  outside  the  context  of  the  Maximum  Entropy  Principle.  The 
more  rigorous  definition  we  adopt  here  is  in  the  form  of  the  modified  Legendre  transform 


g(T,p)  =  min  [e  +  pv  —  T s(e,v)\  .  (18) 

e,  v 

The  definition  (18)  reflects  the  Maximum  Entropy  Principle  applied  to  a  small  thermo¬ 
dynamic  system  in  contact  with  a  reservoir  of  heat  and  work,  and  clearly  demonstrates 
that  the  absolute  temperature  T  and  pressure  p  are  the  variables  on  which  the  Gibbs  free 
energy  naturally  depends. 

The  Gibbs  free  energy  is  also  called  the  free  enthalpy ,  and  the  term  free  energy  is  usually 
reserved  for  the  (specific)  Helmholtz  free  energy  defined  by  the  formula 


f(T,  v )  =  min  [e  —  T  s(e,  u)]  . 

e 


(19) 


The  Helmholtz  free  energy  is  useful  in  describing  a  thermodynamic  system  of  fixed  volume 
(isochoric  process)  attached  to  a  reservoir  of  heat  at  given  temperature  T. 

The  most  important  consequence  of  the  definition  (18)  is  the  fact  that  g(T,p)  is  a  continu¬ 
ous  concave  function  of  T  and  p.  This  conclusion  directly  follows  from  the  observation  that 
a  function  defined  as  the  minimum  of  a  family  of  linear  functions  is  concave.  Continuity 
follows  from  the  general  fact  that  the  minimum  of  a  family  of  continuous  functions  (linear 
functions  in  our  case)  is  a  continuous  function.  Likewise,  the  definition  (19)  implies  that 
f(T,  v )  is  a  continuous  concave  function  of  T  at  any  fixed  v.  To  illustrate  this  just  draw 
a  collection  of  straight  lines  on  a  plane  and  fit  a  curve  below  them,  tangent  to  each  line. 

Assuming  that  g(T,p)  is  a  smooth  function,  we  immediately  arrive  at  the  dependence  of 
the  primary  variables  on  the  potential  variables 


rp  dg  dg 
e  =  9-TBf-pdp' 


dg  =dg_ 
dT  ’  dp' 


(20) 
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Of  course,  all  the  primary  variables  (including  entropy)  in  (20)  are  now  functions  of  T 
and  p.  Note  that,  with  the  help  of  the  Legendre  transform,  the  solution  of  the  system 
of  nonlinear  equations  (14)  is  readily  obtained  by  simple  differentiation.  This  procedure 
automatically  singles  out  the  unique  solution  to  (14)  consistent  with  the  Maximum  Entropy 
Principle.  It  is  worthwhile  noting  that  the  map 


(e,u)  H- 


/ ds  ds\ 
\  de  ’  dv  ) 


(21) 


is  also  called  a  Legendre  transform  (or  a  contact  transformation).  However,  in  this  weak 
form  the  Legendre  transform  is  not  very  useful  in  thermodynamics. 

The  map  (21)  may  not  be  invertible  a  priori.  However,  it  is  invertible  if  the  entropy 
dependence  s  =  s(e,v)  is  a  strictly  concave  function.  The  later  means  that  the  following 
strict  inequalities  take  place 

d2s  d2s  d2s  d2s  /  d2s  \2 

de2  <  ’  dv2  <  ’  de2  dv 2  >  \<9e<9t>  J 


In  this  case  the  Gibbs  free  energy  is  also  a  strictly  concave  function 

2 


d2g  d2g  _ 

— —  <  0  — -  <  0  — —  — -  >  , 

dT2  ’  dp2  ’  dT2  dp2  \dTdp)  ' 


d2g  d2g  f  d2g  V 


(23) 


Moreover,  the  original  entropy  function  s(e,  v )  coincides  with  its  double  Legendre  trans¬ 
form  s(e,  v)  given  by  the  formula 


s(e,  v )  =  min 

T,p 


e  +  pv  -  g(T,p) 


(24) 


In  other  words,  the  Legendre  transform  on  concave  functions  is  involutive:  s(e,  v )  =  s(e,  v). 

From  (23)  a  bunch  of  important  consequences  will  be  now  deduced.  Let  us  introduce  the 
specific  heat  capacity  at  constant  pressure  cp,  isothermal  compressibility  ot,  and  isobaric 
thermal  expansivity  bp  by  the  formulae 

ds  1  dv  1  dv 

Cp  dT  aT  v  dp'  p  v  dT 

These  thermodynamic  properties  are  measurable  in  experiments  aimed  at  determining  the 
equation  of  state  of  a  particular  material.  The  heat  capacity  cp  and  isothermal  compress¬ 
ibility  ot  are  always  positive  since  the  function  g(T,p )  is  strictly  concave 

d2a  1  d2 a 

cr  =  ~T  gp  >  0 .  (26) 

In  contrast  to  this,  nothing  prevents  the  isobaric  thermal  expansivity  bp  to  have  any  sign. 
For  example,  water  at  standard  pressure  condition  (p  =  0.1  MPa)  has  a  local  maximum  of 
density  with  respect  of  temperature  near  T  =  277  K  (4°  C).  However,  due  to  the  identity 

b  =  1  92  9 
p  v  dT  dp 
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the  last  condition  (23)  imposes  the  bound 


\bp\  < 


/ 


OT  cp 

vT  ‘ 


(28) 


The  properties  cp,  ax  and  bp  as  functions  of  T  and  p  are  not  independent  but  must  satisfy 
compatibility  conditions  resulting  from  the  fact  that  all  of  them  are  defined  in  terms  of 
the  single  function  g(T,p).  For  example,  the  following  identities  hold 


d(cP/T)  d(bp  v )  =  d{aTv)  d{bpv)  = 
dp  dT  ~  ’  8T  Bp  ~ 


This  fact  is  very  useful  as  significantly  reduces  the  number  of  experiments  for  determining 
the  equation  of  state  of  a  material  in  question. 

Let  cv  be  the  specific  heat  capacity  at  constant  volume.  It  is  straightforward  to  prove  the 
equality  [Reif  1965] 


Cp  cv 


b2pTv 

(ix 


(30) 


which  particularly  ensures  the  inequality  cp>  cv.  The  strict  concavity  of  g(T,  p)  guarantees 
that  cv  is  positive.  Indeed,  direct  calculations  result  in  the  inequality 


T 

d2g  d2g 

(  d2g  \2' 

ax  v 

dT 2  dp 2 

\  dT  dp  ) 

(31) 


(This  result  can  be  also  obtained  from  the  strict  concavity  of  the  Helmholtz  free  energy 
with  respect  of  T  because 


s 


df_ 

dT  ’ 


Cy 


ds  _  d2f 

Tdf  =  ~Tdr2 


>  0. 


(32) 


Here  the  independent  variables  are  T  and  v.) 

The  positivity  of  ax,  cp  and  cv  is  an  important  consequence  of  the  definition  of  the  Gibbs 
free  energy  in  terms  of  the  Legendre  transform.  In  many  circumstances  this  allows  one  to 
prove  Le  Chatelier’s  principle  which  states  [Reif  1965]: 


If  a  system  is  in  stable  equilibrium,  then  any  spontaneous  change  of  its  parame¬ 
ters  must  bring  about  processes  which  tend  to  restore  the  system  to  equilibrium. 


A  thermodynamically  consistent  mathematical  model  for  real-world  physical  phenomena 
must  respect  this  important  principle  in  the  first  place. 


2.2  Phase  transitions 

Phase  transitions  of  the  first  kind  are  characterized  by  an  abrupt  dependence  of  the  pri¬ 
mary  variables  (e,  v)  on  the  potential  variables  (T,  p) .  If  the  entropy  dependence  s  =  s(e,  v) 
is  a  strictly  concave  function  of  e  and  v,  then  the  Legendre  transform  is  also  a  strictly 
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concave  function,  and  in  particular  there  will  be  an  invertible  mapping  between  ( T,p )  and 
(e,  v )  in  terms  of  the  Gibbs  free  energy. 

It  is  conceivable  from  the  above  that,  within  this  framework,  phase  transitions  can  be 
modelled  under  the  assumption  that  s(e,  v )  is  not  a  concave  (or  strictly  concave)  function 
of  the  primary  variables.  In  this  case  the  function  g  (T,  p )  will  have  discontinuous  partial 
derivatives,  and  therefore  a  discontinuous  dependence  of  the  primary  variables  on  the 
potential  variables  will  be  observed.  The  jumps  of  the  primary  variables  is  nothing  else 
than  the  latent  energy  and  latent  volume  of  phase  transition.  This  kind  of  singularities 
of  Legendre  transforms  is  analysed  in  the  theory  of  catastrophes  [Broker  1975,  Arnold, 
Gusein-Zade  &  Varchenko  1985]. 


Figure  3:  Entropy  versus  energy 


Figure  4:  Free  energy  versus  temperature 

Let  us  illustrate  the  described  discontinuous  behaviour  using  an  example  of  an  incompress¬ 
ible  (one-parametric)  phase  transition.  We  suppress  the  dependence  on  v  and  consider  a 
non-concave  dependence  of  entropy  on  energy  as  shown  in  Figure  3.  The  graph  contains  a 
local  region  of  non-concavity  (the  ‘dimple’).  The  calculated  Helmholtz  free  energy  is  shown 
in  Figure  4.  Obviously,  some  interval  of  energy  e  is  Legendre  transformed  to  a  single  point, 
the  phase-transition  temperature  T*,  where  the  free  energy  has  a  kink.  Subsequently,  en¬ 
ergy  as  a  function  of  temperature  exhibits  a  jump  at  the  phase-transition  temperature  as 
shown  in  Figure  5.  This  jump  is  the  latent  energy.  The  double  Legendre  transform  of 
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Figure  5:  Energy  versus  temperature 


Figure  6:  Double  Legendre  transformed  entropy  versus  energy 


entropy  is  shown  in  Figure  6.  This  function  of  energy  is  concave  but  not  strictly  concave 
as  the  graph  has  a  region  with  constant  slope  which  removes  the  ‘dimple’  of  the  original 
function.  Outside  of  this  region  the  original  and  ‘corrected’  entropy  functions  coincide. 

It  is  worthwhile  emphasizing  that  the  Legendre  transform  maps  the  one-dimensional  inter¬ 
val  in  the  space  of  energies  to  the  zero-dimensional  point  (T*)  in  the  space  of  temperatures. 
In  other  words,  one  dimension  is  lost  under  the  Legendre  transformation  of  a  non-concave 
entropy. 

Likewise,  phase  transitions  in  a  two-parametric  system  kept  in  contact  with  a  reservoir  at 
given  temperature  T  and  pressure  p  occur  at  a  phase-transition  line  in  the  (T,  p)-plane. 
In  other  words,  some  two-dimensional  region  in  the  (e,  u)-plane  is  Legendre  transformed 
to  a  one-dimensional  line  in  the  (T,p)- plane.3  In  practice  this  phase-transition  line  is 
parametrized  either  by  temperature,  p  =  p*(T),  or  by  pressure,  T  =  T*(p).  On  the  vapouri¬ 
sation  line  the  function  p  =  p*(T)  is  called  the  vapour  pressure,  and  on  the  fusion  line  the 
function  T  =  T*(p)  is  the  pressure-dependent  fusion  temperature.  Both  approaches  can  be 
unified  by  assuming  that  the  points  of  the  phase-transition  line  (T*,p*)  are  parametrized 
by  some  convenient  parameter  or  given  implicitly  by  some  equation  p  (T,  p)  =  0  in  the 
(T,  p)-plane.  Of  course,  in  the  case  of  a  vapourisation  line,  this  equation  makes  sense  up 

3 Theoretically,  this  line  can  collapse  to  a  point  in  a  degenerate  case. 
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to  the  critical  state. 

Due  to  the  global  continuity  of  g(T,p),  the  second  law  of  thermodynamics  (15)  is  valid  in 
finite  variations  when  the  phase-transition  line  is  crossed  at  a  point  (T*,p*),  namely 

T*  bs  =  be  +  p*  bv  .  (33) 


Here  5s,  be  and  bv  are  the  latent  entropy,  energy  and  volume,  respectively.  Note  that,  since 
a  phase  transition  is  a  reversible  thermodynamic  process,  the  latent  entropy,  energy  and 
volume  change  their  signs  if  the  phase-transition  line  is  crossed  in  the  opposite  direction. 
To  avoid  this  ambiguity  the  latent  properties  are  uniquely  defined  by  requiring  that  the 
phase-transition  line  is  crossed  in  a  direction  from  solid  to  liquid  or  from  liquid  to  gas 
(vapour)  phases.  The  term  T*  bs  is  the  latent  heat,  and  it  is  natural  to  call  the  product  p*  bv 
the  latent  work. 

Since  the  differential  d g(T,p)  evaluated  at  a  tangent  element  (dT*,dp*)  to  the  phase- 
transition  line  is  the  same  for  both  phases,  the  Clausius-Clapeyron  equation  [Reif  1965] 
is  immediately  derived 

dp*  _  bs 
d  T*  bv 


Equations  (33)  and  (34)  allow  one  to  express  the  latent  energy  and  latent  volume  in  terms 
of  (T*,p*)  and  5s,  namely 


=  (35) 

The  phase-transition  diagram  and  the  latent  heat  are  measured  in  experiments  for  a  range 
of  temperatures,  and  therefore  be  and  bv  are  completely  determined  by  (35).  Of  course, 
since  the  primary  variables  are  discontinuous,  the  properties  cp,  ax  and  bp  exhibit  delta- 
function  singularities  at  the  phase-transition  line. 

When  modelling  fluid  flows  with  phase  transitions,  it  is  more  consistent  to  use  the  double 
Legendre  transform  (24)  beforehand  as  the  fundamental  equation  of  state.  With  entropy 
thus  defined,  it  is  straightforward  to  check  that  the  Jacobian  of  the  transformation  (e,  v)  H y 
( T, p )  is  singular  at  phase  coexistence  states,  i.e. 


det 


d(T,p) 
d(e,  v) 


=  0. 


Indeed,  we  have  the  conditions  for  (non-strict)  concavity 

d2s  ^  d2s  ^  d2s  d2s  /  d2s  \2 
de2  ~  ’  dv 2  ~  ’  de2  dv 2  ~  \dedv) 


(36) 


(37) 


which  tell  us  that  the  Jacobian  (36)  vanishes  if  either  of  inequalities  becomes  an  equality. 

The  images  of  the  points  at  which  the  Jacobian  of  a  transformation  is  singular  are  called 
critical  values.  According  to  Sard’s  theorem  [Sard  1942,  Sternberg  1964]  the  set  of  the 
critical  values  of  any  (differentiable)  map  has  zero  measure.  In  our  case  of  a  two-paranretric 
thermodynamic  system,  the  critical  values  of  the  Legendre  transform  constitute  the  vapour 
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pressure  curve  whose  area  (Lebesgue  measure)  is  obviously  zero.  Sard’s  theorem  provides 
a  mathematical  foundation  for  Gibbs’  phase  rule  [Gibbs  1948]  though  does  not  prove  it 
in  full.  Strictly  speaking,  Gibbs’  phase  rule  is  valid  for  the  general  configuration  only. 
For  example,  it  is  straightforward  to  design  a  degenerate  entropy  dependence  s  =  s(e,  v ) 
which  maps  the  phase  coexistence  domain  in  the  primary  variables  to  a  single  point  in 
the  potential  variables  by  means  of  the  Legendre  transform.  However,  such  an  equation 
of  state  may  model  no  real  fluid. 

The  described  approach  is  the  most  straightforward  and  thermodynamically  consistent 
way  to  model  phase  transitions.  Phase  transitions  discussed  here  so  far  belong  to  the 
class  of  first-order  phase  transitions.  Second-order  phase  transitions  are  characterized 
by  a  continuous  behaviour  of  the  primary  variables  on  the  potential  variables  but  some 
other  derived  properties,  such  as  heat  capacity,  diverge.  In  this  case  there  is  no  latent 
heat  anymore.  Technically,  this  situation  can  be  modelled  by  assuming  that  the  entropy 
s  =  s(e,  v )  is  not  strictly  concave  at  some  set  in  the  (e,  u)-plane  with  empty  interior,  such 
as  a  line  or  a  point. 

Second-order  phase  transitions  also  take  place  in  conventional  fluids  when  pressure  and 
temperature  are  carefully  controlled  in  such  a  manner  that  the  point  (T,  p)  of  the  phase 
diagram  passes  through  the  critical  point  (Tc,pc),  being  transversal  (not  tangent)  to  the 
phase-transition  line.  The  theoretical  background  of  critical  phenomena  is  presented  in 
many  monographs  [Binney  et  al.  2002]. 


2.3  The  van  der  Waals  equation  of  state 


Consider  the  thermal  equation  of  state  of  the  van  der  Waals  fluid 


p(T,  v ) 


RT  a 

v  —  b  v 2 


(38) 


where  R  is  the  specific  gas  constant,  and  a  and  b  are  positive  parameters.  The  specific 
volume  v  must  be  greater  than  b.  Formula  (38)  reduces  to  the  thermal  equation  of  state 
of  an  ideal  gas  when  both  a  and  b  vanish. 

The  critical  state  of  the  van  der  Waals  fluid  is  obtained  from  the  following  equations  to 
be  solved  simultaneously 

_ ,_v)  =  n  d2p(T,  v ) 

dv  ’  dv2 


which  have  the  solution 


Tc  = 


8  a 

27b  R  ’ 


Pc  27 b2  ' 


(40) 


The  above  formulae  are  expressed  for  a  and  b  in  terms  of  the  critical  temperature  Tc  and 
pressure  pc  as  follows 


27  R2T2  ,  RTC 

- -  ,  b  = - - 

64 pc  8  pc 


(41) 
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The  critical  parameters  Tc  and  pc  are  quite  reliably  measured  for  various  materials,  and 
thus  the  van  der  Waals  equation  of  state  is  completely  specified  given  the  specific  gas 
constant  R. 

The  second  law  of  thermodynamics  (15)  requires  that  the  specific  heat  at  constant  volume 
is  a  function  of  temperature  only:  cv  =  cv(T).  This  is  proved  in  exactly  the  same  way  as 
for  an  ideal  gas.  Everything  that  matters  is  the  linearity  of  p(T,  v )  with  respect  of  T.  The 
dependence  cv  =  cv{T)  is  called  the  caloric  equation  of  state. 

The  thermal  and  caloric  equations  of  state  define  entropy  s  =  s(e,v)  up  to  a  constant. 
Indeed,  let  us  calculate  the  anti-derivative 

u(T)  =  J  cv(T)  AT  ,  (42) 


and  then  calculate  its  inverse  function  u  1(e).  The  inverse  function  exists  because  cv  >  0. 
Then  calculate  the  anti-derivative 

(43) 

in  terms  of  which  the  dependence  of  entropy  on  the  primary  variables  is  established 


w(e)  = 


de 


u 


-l 


(e) 


s(e,v) 


sQ  +  R  log 


v  —  b 
va-b 


+  w 


—  w 


(44) 


where  e0,  vQ  and  sa  are  some  reference  parameters.  For  constant  cv  (perfect  fluid),  it 
reduces  to  the  analytic  expression 


s(e,v) 


s0  +  R  log 


v  —  b 
vQ-b 


.  e  +  a/v 

+  cv  log - - — j—  . 

c0  +  a  v0 


(45) 


This  function  is  not  concave  with  respect  of  (e,  v)  as  can  be  checked  by  direct  differentia¬ 
tion,  though  it  is  strictly  concave  with  respect  of  e  for  any  fixed  v. 

In  order  to  compute  the  Gibbs  free  energy  g(T,p),  it  suffices  to  solve  the  following  cubic 
equation 


pv 3  —  (RT  +  bp)  v2  +  av  —  ab  =  0 


(46) 


for  volume  v  at  given  T  and  p,  and  then  find  the  minimum  value  of  the  expression 


cvT - \-pv  -T 

v 


,  D  !  v  ~h  1  T 

s0  +  R  log - -  +  Cv  log 


vQ-b 


eQ  +  a/vQ_ 


(47) 


over  the  real  roots.  This  minimum  value  is  equal  to  g(T,p). 

Let  us  specify  the  van  der  Waals  equation  of  state  for  water.  In  this  case  we  choose 
R  =  462  J/K/kg,  cv  =  1600  J/K/kg,  Tc  =  647  K  and  pc  =  22  MPa.  Strictly  speaking,  cv  is 
quite  different  for  water  and  its  steam  [Lide  1999].  Nevertheless,  to  capture  the  qualitative 
behaviour  of  phase  transitions  in  water,  we  assume  that  it  is  constant  (perfect  fluid). 

Figure  7  displays  the  van  der  Waals  equation  of  state  (38)  for  two  subcritical  and  one 
supercritical  temperature  values.  On  the  supercritical  isotherm,  pressure  monotonically 
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Figure  1:  The  van  der  Waals  equation  of  state 


Figure  8:  Energy-pressure  isotherms 


decreases  with  increasing  specific  volume,  whereas  the  subcritical  isotherms  are  not  mono¬ 
tone.  Moreover,  pressure  becomes  even  negative  that  does  not  make  sense  from  the  phys¬ 
ical  point  of  view.  This  shortcoming  is  corrected  by  replacing  the  analytic  expression  for 
entropy  with  its  double  Legendre  transform  (24).  Alternatively,  we  can  employ  the  Gibbs 
free  energy  and  calculate  the  primary  variables  in  terms  of  the  potential  variables. 

Figure  8  shows  energy-pressure  isotherms.  As  is  expected,  specific  energy  changes  abruptly 
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at  some  pressure  and  given  subcritical  temperature  defining  the  isotherm.  The  latent 
energy  diminishes  with  increasing  temperature  until  it  vanishes  beyond  the  critical  point. 


Figure  9:  Volume-pressure  isotherms 

Figure  9  shows  volume-pressure  isotherms.  In  fact,  after  swapping  the  coordinate  axes, 
these  are  the  pressure-volume  isotherms  (the  van  der  Waals  equation  of  state)  shown  in 
Figure  7  but  obtained  from  the  Gibbs  free  energy  or,  equivalently,  from  the  double  Leg¬ 
endre  transform  of  entropy.  In  contrast  to  the  van  der  Waals  equation  of  state,  these 
isotherms  are  monotonic,  exhibiting  an  abrupt  change  in  specific  volume.  It  is  straight¬ 
forward  to  verify  that  these  isotherms  can  be  obtained  by  applying  Maxwell’s  equal  area 
rule  [Maxwell  1875,  Reif  1965]. 

The  specific  heat  capacities  as  functions  of  temperature  at  standard  pressure  are  shown  in 
Figure  10.  The  delta-function  singularity  of  the  specific  heat  capacity  at  constant  pressure 
(cp)  is  depicted.  Since  we  model  water  and  its  vapour  by  a  perfect  fluid,  the  specific  heat 
capacity  at  constant  volume  (cv)  is  constant.  In  real  fluids  both  heat  capacities  exhibit 
similar  types  of  singularities  which  are  not  captured  by  this  model.  It  is  also  seen  that  the 
vapourisation  temperature  is  greatly  underpredicted.  This  is  obviously  a  shortcoming  of 
the  van  der  Waals  equation  of  state. 

The  vapour  pressure  curve  for  the  van  der  Waals  fluid  is  shown  in  Figure  11.  This  is 
the  phase  diagram  in  potential  variables.  As  is  expected,  the  vapourisation  line  ends  at 
the  critical  point  ( TClpc )  depicted  by  the  filled  circle.  The  phase  diagram  in  the  plane  of 
the  primary  variables  is  shown  in  Figure  12.  The  coexistent  states  of  the  liquid- vapour 
phases  correspond  to  the  shaded  region.  The  slopes  of  the  straight  lines  are  given  by 
the  relations  (35).  This  two-dimensional  region  of  coexistent  phases  is  mapped  to  the 
one-dimensional  vapourisation  line  by  the  Legendre  transform.  The  critical  state  in  the 
primary  variables  mapped  to  ( Tc,pc )  by  the  Legendre  transform  is  depicted  by  the  filled 
circle. 
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Figure  10:  Specific  heat  capacities 


Figure  11:  Vapour  pressure  curve 


The  latent  heat,  energy  and  work  as  functions  of  the  absolute  temperature  are  shown  in 
Figure  13.  Of  course,  the  sum  of  the  latent  energy  and  work  is  equal  to  the  latent  heat  as 
dictated  by  (33).  All  these  properties  vanish  at  the  critical  temperature  Tc  and  beyond. 

The  van  der  Waals  equation  of  state  does  not  predict  well  the  phase-transition  temperature 
and  pressure  of  a  real  fluid  from  the  known  specific  gas  constant  R  and  critical  state  (: Tc,pc ), 
so  its  application  to  the  modelling  of  real-world  phenomena  is  limited.  Nevertheless,  this 
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Figure  12:  Phase  diagram  in  primary  variables 


Figure  13:  Latent  heat,  energy  and  work 


equation  of  state  is  an  excellent  model  for  theoretical  studies,  and  provides  a  reasonable 
model  for  near  critical  behaviour  of  pure  liquids  and  their  vapours. 
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3  Mathematical  models 

The  modelling  of  phase  transitions  from  the  first  principles  is  based  on  Thermodynamics 
and  Continuum  Mechanics.  In  this  section  we  formulate  the  most  general  model,  and 
discuss  several  approximations. 


3.1  Conservation  laws 

The  conservation  laws  of  continuum  mechanics  constituting  balance  of  mass,  momentum 
and  energy  are  written  in  the  form  [Truesdell  1977] 

||  +  V-M  =  0,  (48) 


d  (pv) 
dt 


+  V  •  (p  v  <g>  v  +  P) 


(49) 


C/o  — * 

—  +V-(ev  +  P-v  +  q)  =  pf-v.  (50) 

Here  t  is  time,  v  velocity,  p  mass  density,  e  energy  density,  P  momentum  flux,  q  heat  flux, 
and  /  external  force.  The  symbol  <8>  denotes  the  tensor  product,  the  centred  dot  stands 
for  the  scalar  (inner)  product,  and  V  is  the  gradient  operator.  The  volumetric  density  of 
energy  is  the  sum  of  the  internal  and  kinetic  energies,  namely 

£  =  p(e+\\v  I2)  (51) 


where  e  is  the  specific  internal  energy.  The  balance  of  angular  momentum 


d  (x  x  p  v) 
dt 


+  V  ■  (x  x  pv  <S>  v  +  x  x  P) 


xx  pf 


(52) 


requires  the  momentum  flux  P  to  be  a  symmetric  second-rank  tensor.  Here  x  is  the 
position  vector  and  the  symbol  x  denotes  the  vector  product.  Note  that  the  negative  of  P 
is  the  stress  tensor. 


In  order  to  complete  the  governing  equations,  we  need  to  resort  to  thermodynamics  and 
rheology.  Thermodynamics  of  compressible  fluids  unambiguously  defines  the  absolute 
temperature  T  and  pressure  p.  For  fluids  at  equilibrium  the  momentum  flux  P  is  expressed 
in  terms  of  pressure  p  by  the  formula 


P  =  pl  (53) 

where  I  is  the  identity  tensor.  In  other  words,  P  ■  n  =  ph  for  any  unit  vector  n,  which  is 
the  manifestation  of  Pascal’s  law.  Of  course,  at  equilibrium  the  heat  flux  vanishes:  q  =  0. 

Rheological  constitutive  equations  are  applicable  to  irreversible  processes  due  to  heat  flux 
and  fluid  motion.  The  momentum  flux  P  and  heat  flux  q  are  determined  in  terms  of  the 
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gradients  of  velocity  and  temperature.  Neglecting  reciprocal  cross-effects,  the  momentum 
flux  depends  on  the  deformation-rate  tensor  D  =  (Vu)sym  according  to  Stokes’  law 

P  =  pi  —  2pD  (54) 

where  p  is  the  dynamic  coefficient  of  viscosity  (p  is  also  called  the  shear  coefficient  of 
viscosity  to  distinguish  it  from  the  dilatational  coefficient  which  is  usually  neglected  espe¬ 
cially  in  nearly  incompressible  flow).  The  heat  flux  depends  on  the  temperature  gradient 
according  to  Fourier’s  law 

q  =  -n\7T  (55) 


where  n  is  thermal  conductivity.  The  physical  properties  p  and  n  depend  on  the  thermo¬ 
dynamic  state  of  material. 

We  assume  the  material  is  determined  by  a  given  dependence  of  specific  entropy  s  on 
the  specific  internal  energy  e  and  specific  volume  v  =  1/p.  This  assumption  is  in  fact 
the  definition  of  ffuids.  A  sufficiently  small  control  volume  of  fluid  with  nearly  uniform 
distribution  of  mass  and  energy  is  identified  with  a  thermodynamic  system  attached  to  the 
environment,  a  reservoir  of  heat  and  work  with  the  absolute  temperature  T  and  pressure  p 
evaluated  at  the  location  of  the  control  volume. 

We  can  express  the  absolute  temperature  and  pressure  in  terms  of  the  double  Legendre 
transformed  entropy 


T(e,v)  =  1 


ds(e ,  v ) 
de 


p(e,v) 


ds(e,v)  /ds(e,v) 
dv  /  de 


(56) 


It  is  important  to  emphasize  that  the  original  function  s  =  s(e,v)  is  not  suitable  for 
determining  the  absolute  temperature  and  pressure,  unless  it  is  strictly  concave  everywhere 
and  hence  coincides  with  s(e,v).  Alternatively,  we  can  express  the  primary  variables  in 
terms  of  the  Gibbs  free  energy 


e(T,p)  =  g(T,p)  —  T 


dg(T,p) 

dT 


~P 


dg(T,P) 

dp 


p(T,  p)  =  1 


j  dg(T,  p) 
/  dp 


(57) 


The  choice  between  the  representations  (56)  and  (57)  for  completing  the  mathematical 
problem  depends  on  whether  one  needs  to  use  the  primary  or  potential  variables  as  inde¬ 
pendent  variables. 


3.2  Conditions  at  a  strong  discontinuity 

The  conservation  laws  (48)-(50)  are  applicable  to  any  flow  of  a  continuum  medium,  even 
to  flows  containing  strong  discontinuities  such  as  shocks  and  phase-change  fronts.  In  this 
case  the  divergence  operator  has  to  be  treated  in  the  generalized  sense. 

The  conservation  laws  of  mass,  momentum  and  energy  result  in  the  following  conditions 
at  the  strong  discontinuity  assumed  to  be  a  smooth  surface 

\p(v-n-Vn)\  1  =  0,  (58) 
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[pv  (v-  n  -  Vn)  +  P  ■  n]t  =  0  , 


(59) 


[e  (v  ■  n  —  Vn)  +  v  ■  P  ■  n  +  q  •  n]l  =  0  ,  (60) 

where  n  is  a  unit  normal  vector  to  the  surface,  Vn  is  the  speed  of  propagation  of  the  surface 
along  the  normal  vector,  and  the  brackets  [•]  +  denote  the  jump  operation 


[f]t  =  /+-/-,  f±{x,t)=  lim  f(x  +  zn(x,t),t). 


(61) 


For  example,  if  the  moving  strong  discontinuity  is  described  implicitly  by  the  equation 
F(x,  t)  =  0,  then 


VF  _  dF/dt 

~  |VF|  ’  n  “  ~  |VF| 


(62) 


By  the  way,  the  latter  formulae  provide  the  continuation  of  n  and  Vn  to  the  vicinity  of  the 
strong  discontinuity,  so  the  slightly  abused  notation  for  the  jump  operation  in  (58)— (60) , 
involving  n  and  Vn  defined  at  the  surface  only,  makes  sense. 

As  particular  cases  the  conditions  (58)-(60)  produce  the  Rankine-Hugoniot  conditions  at 
a  shock  front  [Courant  &;  Friedrichs  1948]  and  the  Stefan  condition  at  a  phase-change 
front  [Meirmanov  1992],  The  Stefan  condition  is  derived  by  assuming  zero  velocity.  This 
immediately  demands  that  density  is  continuous  across  the  phase-change  front:  [p\+  =  0. 
Therefore,  since  p  =  p+  =  p_,  the  energy  balance  condition  (60)  reduces  to  the  Stefan 
condition 


[e]±  Vn  =  p  [e]+  Un  =  [q-n]t  .  (63) 

In  the  classical  formulation  of  the  Stefan  problem  that  involves  tracking  the  phase-change 
front,  the  absolute  value  of  the  jump  of  the  specific  internal  energy  e  across  the  phase- 
change  front  is  equal  to  5e,  the  specific  latent  heat  of  phase  transition. 


3.3  The  Stefan  problem 

The  general  governing  equations  reduce  to  the  Stefan  problem  by  assuming  zero  velocity 
and  constant  density.  Such  a  situation  arises  when  the  solid  phase  is  at  rest  and  an 
incompressible  liquid  phase  is  of  the  same  density  as  the  solid  phase  and  does  not  move 
during  the  process  of  solidification  or  fusion. 

The  mass  and  momentum  conservation  laws  are  automatically  satisfied,  so  the  only  non¬ 
trivial  equation  is  the  balance  of  energy  which  reduces  to 

^  =  V  ■  (kVT)  (64) 
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where  e  =  p  e  with  constant  density  p.  The  Stefan  condition  (63)  along  with  the  condition 
for  fusion  temperature  is  specialized  to  the  equations  at  the  phase-change  front 

+ 

=  0  ,  T+  =  T_  =  T*  .  (65) 


[e]t  Vn  + 


dT 

dn 


Assuming  given  dependences  e  =  e(T)  and  n  =  k(T),  we  end  up  with  the  classical  for¬ 
mulation  of  the  Stefan  problem  in  which  the  location  of  the  phase-change  front  is  part  of 
solution. 

There  are  many  numerical  schemes  aimed  at  an  accurate  solution  to  the  classical  Stefan 
problem,  however  their  application  to  real-world  problems  is  limited.  This  is  because  the 
tracking  of  the  phase-change  front  is  tedious  as  the  front  can  change  topology  with  evolving 
time.  Indeed,  the  fronts  can  grow  from  nothing  (new  phase  nucleation),  merge,  break  and 
disappear  altogether.  So,  it  is  desirable  to  provide  a  uniform  formulation  of  the  Stefan 
problem  which  does  not  involve  tracking  the  phase-change  fronts  explicitly,  and  can  be 
employed  in  Computational  Fluid  Dynamics  (CFD)  simulation  of  real-world  problems  in 
a  straightforward  manner. 


3.3.1  Fictitious  capacity  method 

Let  us  introduce  the  volumetric  heat  capacity  p  =  pcp  which  is  equal  to  the  derivative  of 
the  volumetric  internal  energy 

MT) 

dT 

The  energy-balance  equation  is  specialized  to  the  classical  heat  conduction  equation  with 
variable  coefficients 

p(T)^=V-[k(T)VT)}.  (67) 

When  phase  transitions  take  place  the  volumetric  heat  capacity  exhibits  a  delta-function 
singularity  at  T  =  T*. 

In  order  to  avoid  this  problem,  we  approximate  the  coefficient  p(T)  in  (67)  by  the  centred 
finite  difference 

e  (T  +  6T/2)  -  e  (T  -  ST/2) 

5T 

where  5T  is  some  regularization  parameter.  It  is  natural  to  call  this  regularization  the 
fictitious-capacity  approximation.  The  advantage  of  this  formulation  is  in  the  fact  that 
there  is  no  need  to  track  the  phase-transition  front.  It  is  conceivable  that  the  regularization 
parameter  6T  could  be  neither  very  small  nor  very  large,  so  there  arises  the  problem  of 
an  optimal  choice  of  5T. 

Many  efficient  numerical  algorithms  are  designed  for  solution  of  the  problem  (67),  however, 
all  of  them  are  not  conservative  in  nature  because  the  balance  of  energy  is  not  approx¬ 
imated  exactly.  Due  to  the  fictitious-capacity  regularization,  implicit  numerical  schemes 
can  be  applied  to  this  class  of  problems.  Finite  Element  Method  (FEM)  is  a  good  choice 
for  numerical  solution  of  the  Stefan  problem  in  temperature  formulation. 


UNCLASSIFIED 


21 


DST-Group-TR-3139 


UNCLASSIFIED 


3.3.2  Enthalpy  formulation 

The  alternative  ‘enthalpy’  formulation  is  in  terms  of  energy 

^  =  V-Ke)VT(e)]  (69) 

where  the  thermal  conductivity  n(e)  is  now  a  given  function  of  energy.  Note  that  this 
formulation  requires  a  suitable  interpolation  of  the  thermal  conductivity  k  over  the  range 
of  energies  corresponding  to  the  thermodynamic  states  of  phase  coexistence,  namely,  where 
T(e)  =  T*.  Of  course,  the  previous  temperature  dependence  k(T )  can  be  superposed 
with  T{e)  to  obtain  the  desired  dependence  but  this  boils  down  to  determining  k  at  the 
fusion  temperature  T*. 

The  enthalpy  formulation  allows  one  to  design  conservative  numerical  schemes.  Usually, 
such  schemes  are  explicit  as  the  problem  degenerates  at  T(e)  =  T*.  Finite  Volume  Method 
(FVM)  is  well  suited  to  this  formulation. 


3.3.3  Exact  similarity  solution 


The  exact  solution  to  the  Stefan  problem  we  are  going  to  describe  has  the  following  physical 
interpretation.  At  time  t  =  0  two  phases  of  a  substance  occupy  two  half-planes,  {x  >  0} 
and  {x  <  0},  are  kept  at  constant  temperatures  T±  .  and  have  temperature-independent 
thermal  conductivities  k±  and  volumetric  heat  capacities  rj±.  Here  the  ‘plus’  and  ‘minus’ 
subscripts  correspond  to  the  positive  and  negative  values  of  x,  respectively.  The  phase- 
transition  temperature  T*  is  implied  to  be  between  T£°  and  T“.  The  latent  heat  of  fusion 
per  unit  volume  is  Se  =  pde. 

The  front  position  is  defined  by  the  law  x  =  2A  yft  with  some  A  to  be  determined.  The 
following  similarity  solution  satisfies  the  heat  conductivity  equation  [Hill  1987] 


T(x,  t)  =  T±  +  (T*  -  T|°) 


erfc  ( ± 


vw«±  c) 


e  = 


erfc  (±\Jti±/k±  A^  2\/t 


(70) 


where  the  ‘plus’  and  ‘minus’  subscripts  are  now  taken  for  x  >  2A  y/t  and  x  <  2A  y/t, 
respectively,  and 


erfc(x) 


OO 

J  exp  (— x 2)  dx 

X 


(71) 


is  the  complementary  error  function.  For  t  >  0,  the  temperature  T(x,t )  is  continuous 
everywhere,  tends  to  T±  as  x  — »  Too,  and  approaches  T*  at  the  phase-change  front.  It 
is  straightforward  to  check  that  the  Stefan  condition  at  the  phase-change  front  x  =  2A  y/t 
results  in  the  transcendent  equation 


sgn  (T“ 


{T+  ~  T '*)  V71+  K+  +  iT-  ~  T *)  y/V-  K-  =  Q 

T  (+y/rj+/K+  A)  $  yJr}-/K-  A^ 


(72) 
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where  <h(x)  =  ypK  exp  (x2)  erfc(x),  the  scaled  complementary  error  function  with  asymp¬ 
totic  behaviour  <J>(x)  «  1/x  as  x  — >  +oo. 

The  transcendent  equation  (72)  has  a  unique  solution  A  which  can  be  obtained  numerically 
in  an  efficient  way.  This  exact  solution  is  used  as  a  benchmark  to  test  numerical  schemes. 


3.3.4  Verification  of  uniform  numerical  schemes 

Verification  is  the  process  of  evaluating  a  numerical  scheme  against  a  benchmark  solution. 
This  process  should  not  be  confused  with  validation  which  is  the  process  of  testing  a 
model  against  real-world  experiments.  A  particular  model  can  pass  verification  but  fail 
validation,  for  example,  because  the  underlying  mathematical  model  may  not  include  all 
the  essential  physical  effects. 


Figure  14-'  Entropy  and  temperature  of  the  ice-water  system  versus  volumetric  energy 

The  exact  similarity  solution  is  used  to  verify  the  explicit  conservative  scheme  and  the  im¬ 
plicit  fictitious-capacity  scheme.  The  material  is  the  ice-water  system  at  standard  pressure 
with  the  following  physical  properties  [Lide  1999]: 

•  fusion  (melting)  temperature  T*  =  273  K, 

•  volumetric  latent  heat  Se  =  334MJ/m3, 

•  volumetric  heat  capacity  of  ice  r]s  =  2.05MJ/m3/K, 

•  volumetric  heat  capacity  of  water  ly  =  4.22MJ/m3/K, 

•  thermal  conductivity  of  ice  ns  =  2.18W/m/K, 

•  thermal  conductivity  of  water  m  =  0.56W/m/K. 
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Time  Is]:  3600 


Figure  15:  Conservative  scheme 


Time  (&]:  3600 


Figure  16:  Fictitious- capacity  scheme:  5T  =  IK 
The  dependences  s  =  s(e)  and  T  =  T(e)  are  shown  in  Figure  14. 

In  numerical  simulations  the  following  values  for  the  parameters  defining  the  benchmark 
solution  are  taken:  K-  =  m,  k+  =  ks,  r/_  =  rji,  rj+  =  rjs,  T^°  =  350 K  and  X+°  =  270 K. 
The  spatial  interval  is  {0  m  <  X  <  0.1  m},  the  phase-transition  front  is  located  at  the 
middle  x  =  0.05  m,  and  the  initial  temperature  is  piecewise  constant  (this  corresponds  to 
t  =  0  of  the  similarity  solution).  The  Dirichlet  boundary  conditions  at  both  ends  of  the 
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Time  (&]:  3600 


Figure  17:  Fictitious- capacity  scheme:  5T  =  2K 


Figure  18:  Fictitious- capacity  scheme:  5T  =  4K 


spatial  interval  are  maintained  according  to  the  values  given  by  the  exact  solution.  The 
spatial  interval  is  discretized  by  100  evenly  spaced  points  for  both  explicit  conservative 
and  implicit  fictitious-capacity  schemes. 

Figure  15  shows  the  numerical  solution  obtained  with  the  conservative  scheme.  It  is  seen 
a  good  approximation  of  the  temperature  field  and  satisfactory  resolution  of  the  position 
of  the  phase-change  front  at  the  final  time. 
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Figures  16-18  show  the  numerical  solution  obtained  with  the  implicit  fictitious-capacity 
scheme  (with  timestep  r  =  Is)  for  ST  =  1  K,  ST  =  2K  and  ST  =  4K,  respectively.  As  a 
general  fact,  non-conservative  schemes  do  not  predict  the  position  of  a  strong  discontinuity 
well,  and  this  is  evident  for  ST  =  1  K  and  particularly  for  ST  =  4K.  The  optimal  choice  of 
the  regularization  parameter  is  ST  =  2  K  but  this  choice  is  not  universal  as  it  may  change 
in  other  configurations. 


3.4  Models  for  two-parametric  phase  transitions 

The  most  general  model  of  phase  transitions  in  fluids  written  in  the  primary  (conservative) 
variables  directly  utilizes  the  conservation  laws  (48)-(50)  where  the  potential  variables 
given  by  (56)  appear  in  the  constitutive  equations  (54)  and  (55).  The  main  difficulty  is 
associated  with  the  fact  that  the  mapping  (56)  degenerates  at  the  phase-coexistence  states, 
and  therefore  only  explicit  algorithms  can  be  effectively  applied.  There  is  also  a  challenge 
to  implement  the  double  Legendre  transform  of  entropy,  and  eventually  produce  (56), 
which  has  to  be  correlated  with  the  experimentally  obtained  vapour  pressure  curve  and 
latent  heat  of  the  material  in  question. 

Finite  Volume  Method  (FVM)  is  well  suited  to  tackle  this  type  of  problem.  The  implemen¬ 
tation  of  an  explicit  conservative  scheme  is  straightforward.  Discretize  the  computational 
domain  into  a  collection  of  control  volumes  (cells)  separated  by  facets  oriented  by  normal 
vectors,  assign  primary  variables  to  each  cell  and  calculate  normal  fluxes  between  neigh¬ 
bouring  cells.  Choose  a  timestep  small  enough  to  preserve  numerical  stability  and  update 
the  solution  at  the  next  time  level.  The  main  problem  is  associated  with  an  accurate 
evaluation  of  fluxes  between  two  cells  in  contact  without  introducing  numerical  instabil¬ 
ity.  Such  conservative  schemes  are  uniform  as  there  is  no  need  to  track  the  phase-change 
fronts. 

A  uniform  scheme  can  be  designed  in  potential  variables  as  well.  However,  such  schemes  are 
not  conservative  in  nature.  Using  (57)  and  the  definitions  of  heat  capacity,  compressibility 
and  expansivity,  the  governing  equations  take  the  form 

°T  (ft  +  W)  ~  bP  VT)  +  V  •  v  =  0,  (73) 

P  (^  +  F-Wj  +VP-V'(2^)=P/,  (74) 


p  cp 


+  V 


VTJ 


Tbp 


+  v  ■  Vp 


V  •  (k  VT)  +  2  jjD  :  D  . 


(75) 


As  in  the  Stefan  problem,  the  main  difficulty  is  caused  by  the  discontinuous  depen¬ 
dences  (57)  which  results  in  delta- function  singularities  of  the  heat  capacity  cp.  com¬ 
pressibility  ot  and  expansivity  bp. 

Following  the  spirit  of  the  fictitious-capacity  method,  we  approximate  these  coefficients 
by  the  finite  differences 


Cp 


s(T  +  ST/2,p)~s{T-ST/2,p) 
ST 


(76) 
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ClT 


log  [v(T,p-  dp/2)/  v  (T,  p  +  dp/ 2)] 
Sp 


(77) 


— 


log \v(T  +  6T/2,p)/v(T  -  5T/2,p)] 
5T 


(78) 


where  and  5p  are  some  regularization  parameters  of  dimension  of  temperature  and 
pressure,  respectively,  and  s(T,p )  and  v(T,p )  are  expressed  in  terms  of  g(T,p).  The 
problem  of  an  optimal  choice  for  the  regularization  parameters  arises. 


3.4.1  Lubrication  theory  approximation 


The  system  of  equations  (73)-(75)  for  the  potential  variables  (T,p)  and  velocity  v  is  of 
a  mixed  type.  In  the  lubrication  theory  approximation  this  system  of  equations  can  be 
reduced  to  a  parabolic  system  of  equations.  This  approximation  is  valid  for  fluid  flow  in 
thin  gaps  and  tubes  as  well  as  in  a  porous  medium. 

For  the  sake  of  simplicity,  let  us  neglect  the  external  force  /.  For  a  slow  viscous  flow,  the 
solution  to  the  momentum  balance  equation  can  be  approximated  by  Darcy’s  law 

v  =  —(J)X7p  (79) 

where  fluidity  4>  satisfies  the  Poisson  equation 


V  •  (ji  V0)  +  1  =  0  . 


(80) 


Equation  (80)  generalizes  the  Hele-Shaw  flow.  Due  to  adherence  condition  for  velocity, 
the  fluidity  </>  must  vanish  at  contact  with  a  rigid  body,  such  as  a  container  wall.  In 
porous  media,  the  fluidity  cj)  becomes  a  given  coefficient  of  filtration  after  the  procedure  of 
homogenization.  In  this  case  the  conservation  laws  must  be  modified  to  include  porosity. 

With  this  simplification  the  conservation  laws  of  energy  and  mass  reduce  to  the  system  of 
convection-diffusion  equations 


{  du  ■  \ 

Cij  (  +  V  ■  Vuj  J  =  V  •  (Kij  Vuj)  +  Fi 

where 


[  d2g 

d2g  1 

r  k 

1 

dT2 

dT  dp 

,  K  = 

T 

0 

d2g 

d2g 

0 

0 

dT  dp 

dp2 

r  t  i 

'  2 pD : D  k  VT  2  ' 

u  = 

,  F  = 

' J ’  1  'J12 

p 

0 

(81) 


(82) 


(83) 
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and  the  velocity  v  is  expressed  by  Darcy’s  law  (79).  The  ‘capacity’  matrix  C  and  ‘conduc¬ 
tivity’  matrix  K  are  symmetric  and  positive  definite,  and  hence  the  governing  equations 
form  a  parabolic  system  of  equations. 

This  system  of  equations  is  not  standard.  Indeed,  the  fluidity  4>  is  part  of  the  solution  and 
must  satisfy  (80)  with  a  zero  boundary  condition  at  the  solid-fluid  interface.  In  particular, 
the  matrix  K  degenerates  near  this  boundary,  causing  problems  in  numerical  simulation. 
Also,  the  matrix  C  has  delta-function  singularities  at  the  phase-transition  curve  and  has 
to  be  regularized,  for  example,  by  centred  differences.  This  procedure  results  in  large 
variation  of  the  capacity  coefficients  that  may  cause  problems  as  well.  Finally,  the  velocity 
field  is  not  given  but  depends  on  the  solution. 

The  advantage  of  this  symmetric  form  of  the  parabolic  system  of  equations  is  in  the  fact 
that  C  and  K  can  be  diagonalized  simultaneously.  This  fact  is  important  for  designing  a 
stable  numerical  scheme  for  diffusion-convection  flow  using  artificial  balancing  diffusion. 
The  algorithm  is  based  on  the  solution  of  the  characteristic  equation 

\C  —  \K\  =  0  (84) 


for  generalized  eigenvalues  A,  and  then  finding  the  matrix  of  generalized  eigenvectors  E. 
It  is  guaranteed  that  we  get  two  positive  generalized  eigenvalues  from  which  we  form  the 
diagonal  matrix 


A 


A+  0 
0  A_ 


and  represent 


C  =  KEAE-1 . 


(85) 


(86) 


The  artificial  balancing  diffusion  matrix  is  given  by  the  formula 
A  =  KES( A)  E~l 


(87) 


where 


<5  (A) 


5 (Pe+)  0 

0  5 (Pe_) 


(88) 


Here  5{x )  =  x  coth(x)  —  1  is  the  transcendent  function  appearing  in  scalar  convection- 
diffusion  equation,  Pe-t  =  h  X±  |u|  are  the  local  Peclet  numbers,  and  h  is  the  finite  element 
size  in  the  direction  of  velocity.  Anisotropic  diffusion  is  controlled  by  the  dyad  v  <8>  v  to 
reduce  excessive  cross-wind  numerical  diffusion.  This  procedure  generalizes  the  Streamline 
Upwind  Petrov-Galerkin  (SUPG)  method  [Donea  &  Huerta  2003]  to  convection-diffusion 
systems  of  equations. 


3.4.2  Simulation  results 

Consider  a  one-dimensional  viscous  flow  of  fluid  modelled  by  the  van  der  Waals  equation 
of  state,  in  the  lubrication  theory  approximation.  The  initial  state  of  the  fluid  in  the  liquid 
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phase  is  uniform  (constant  temperature  and  pressure),  and  the  flow  is  caused  by  a  gradual 
increase  of  temperature  at  the  left  solid  wall.  The  right  wall  is  permeable  and  insulated, 
where  pressure  is  kept  at  the  initial  reference  pressure,  so  the  fluid  will  start  moving  in 
the  right  direction  due  to  a  pressure  gradient. 


0  0.1  0.2  0.3  0.4  0,5  0.8  0,7  08  0.9  1 

Distance 


Figure  19:  Temperature  and  over-pressure. 


Figure  20:  Specific  energy  and  volume. 

We  employ  dimensionless  variables  in  the  van  der  Waals  equation  of  state  which  correspond 
to  R  =  1,  cv  =  10,  Tc  =  1  and  pc  =  1 .  The  initial  reference  temperature  and  pressure 
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Figure  21:  Fluidity  and  velocity. 


are  equal  to  To  =  0.6  and  po  =  0.1,  respectively.  The  viscosity  and  thermal  conductivity 
depend  on  the  specific  volume  according  to  the  formulae  p,(v)  =  0.01  \/vq/v  and  k(v)  = 
0.1  \Jvq/v  where  vq  is  the  specific  volume  at  this  reference  state.  At  the  left  boundary 
x  =  0  of  the  spatial  interval  {0  <  x  <  1}  the  onset  temperature  is  imposed  which  linearly 
increases  from  the  initial  value  To  to  T\  =  0.8  within  the  time  interval  0  <  t  <  1  and  then 
kept  constant  up  to  the  final  time  t  =  100.  At  the  right  boundary  x  =  1,  pressure  is  equal 
to  po  and  the  Neumann  condition  (zero  heat  flux)  is  imposed  for  temperature.  For  fluidity, 
the  Dirichlet  condition  4>  =  0  is  imposed  at  x  =  0  and  the  Neumann  condition  at  x  =  1. 
Note  that  the  Neumann  condition  is  natural  in  FEM.  The  spatial  interval  is  discretized 
by  1000  points,  and  the  regularization  parameters  are  equal  to  5T  =  0.01  and  5p  =  0.01. 


A  fully  vectorized  FEM  code  is  developed  in  MATLAB®.  In  order  to  achieve  numerical 
stability,  second-order  integration  in  time  is  implemented  with  iterations  to  update  the 
wildly  varying  nonlinear  coefficients  of  the  governing  equations.  The  briefly  described 
SUPG  method  is  used  to  add  artificial  balancing  diffusion  to  the  discretized  transport 
terms  in  the  parabolic  system  of  equations  (81). 


The  distribution  of  temperature  and  over-pressure  (p  —  po)  are  shown  in  Figure  19,  and 
the  specific  energy  and  volume  are  depicted  in  Figure  20.  As  is  expected,  the  potential 
variables  are  continuous  whereas  the  primary  variables  exhibit  a  discontinuous  behaviour. 
The  fluidity  and  velocity  shown  in  Figure  21  are  also  continuous,  though  velocity  exhibits 
quite  a  steep  slope  near  the  phase-transition  front. 


30 


UNCLASSIFIED 


UNCLASSIFIED 


DST-Group-TR-3139 


4  Discussion 

In  Continuum  Mechanics  the  primary  variables  appear  in  the  conservation  laws  of  mass 
and  energy,  whereas  entropy  is  a  new  concept  that  leads  to  the  definition  of  the  absolute 
temperature  and  pressure.  In  this  paper  we  apply  a  consistent  methodology  for  mod¬ 
elling  phase  transitions  using  the  Gibbs  free  energy  defined  by  the  Legendre  transform 
of  entropy  (18)  as  opposed  to  the  traditional  definition  (17).  This  approach  immediately 
provides  us  with  a  bunch  of  important  consequences,  such  as  the  global  continuity  and 
concavity  of  the  Gibbs  free  energy  that  resulted  in  the  positivity  of  specific  heat  capacities 
and  compressibility,  and  the  relations  (35)  for  the  latent  energy  and  volume  in  terms  of 
the  latent  entropy  and  vapour  pressure  curve. 

In  many  textbooks,  fluids  are  defined  as  substances  whose  stress  tensor  at  equilibrium 
is  spherical  [Batchelor  1967,  Prandtl  1967],  or,  equivalently,  for  which  Pascal’s  law  (53) 
widely  used  in  hydraulics  is  valid.1  Methodologically,  it  is  more  natural  to  define  fluids 
as  substances  described  by  a  two-parametric  thermodynamic  system  taking  the  specific 
internal  energy  and  volume  as  primary  variables.  All  the  distinguishing  properties  of 
fluids  can  be  deduced  from  the  known  dependence  s  =  s(e,  v )  using  the  Maximum  Entropy 
Principle.  It  is  worthwhile  noting  that  capillary  fluids  can  be  modelled  by  assuming  that 
entropy  additionally  depends  on  the  gradient  of  density  [de  Sobrino  1976,  Antanovskii 
1996]  or  some  other  order  parameter. 

A  priori  the  entropy  dependence  on  the  primary  variables,  s  =  s(e,v),  can  be  quite  arbi¬ 
trary.  This  is  the  fundamental  equation  of  state  defining  a  particular  substance.  However, 
the  dependence  g  =  g(T,p)  must  always  be  a  continuous  concave  function  as  the  mani¬ 
festation  of  the  Maximum  Entropy  Principle,  and  it  will  be  a  serious  mistake  to  model  a 
particular  phenomenon  assuming  this  function  to  be  given  in  an  arbitrary  fashion. 

It  is  emphasized  that,  when  modelling  fluid  flows  with  phase  transitions,  it  is  more  consis¬ 
tent  to  use  the  double  Legendre  transform  (24)  beforehand  as  the  fundamental  equation 
of  state.  This  is  particularly  useful  when  modelling  phase  transitions  numerically  using 
an  explicit  conservative  scheme. 

The  final  remark  is  related  to  the  fact  that  the  ‘capacity’  and  ‘conductivity’  matrices 
are  symmetric  and  positive  definite.  This  important  observation  allows  us  to  symmetrize 
the  governing  equations  in  the  lubrication  theory  approximation  similarly  to  the  sym- 
metrization  procedure  applied  to  the  thermo-diffusive  Stefan  problem  which  models  phase 
transitions  in  materials  with  a  small  concentration  of  impurity  [Antanovskii  1992].  In  this 
case  the  absolute  temperature  and  chemical  potential  are  the  potential  thermodynamic 
variables.  Symmetrization  of  the  governing  equations  of  heat  transfer  relies  on  a  sym¬ 
metric  form  of  the  ‘conductivity’  matrix  which  is  guaranteed  by  the  Minimum  Entropy 
Production  Principle.  Since  two  symmetric  positive  definite  matrices  can  be  simultane¬ 
ously  diagonalized,  many  results  obtained  for  a  scalar  convection-diffusion  equation  are 
extended  to  matrix  convection-diffusion  equations. 


4This  is  true  when  the  effects  of  capillarity  are  neglected. 
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